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Abstract 

The viscous driven cavity problem is solved using a stream-function and vor- 
ticitv formulation for the incompressible Navier-Stokes equations. Analytical 
development and results are presented by Wood[l]. Convergence 1 is accelerated 
by employing grid sequencing and a Y-cycle multigrid relaxation. The formu- 
lation is second-order accurate in space and first-order implicit in time. Tin 1 
coefficients are lagged to loosely couple 1 the 1 stream-function and vorticity equa- 
tions. The linearized system is relaxed with symmetric Gaufi-Seidel. 



‘Aerospace Technologist, Aerothermodvnamics Branch. Aero- and Gas-Dynamics Division, 
w . a. voodQl arc .nasa.gov (757) 801-8355. 
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Introduction 


The present document is the Fortran source* code for solving the incompress- 
ible Navier-Stokes equations inside a two-dimensional driven cavity. The govern- 
ing equations are cast in a streain-function/vorticity formulation as presented 
by Woodfl]. 

Results presented in Ref. [1] focus on streamline patterns, stability issues, 
and multigrid performance. Additional results presented here provide numerical 
tabulations and comparisons with results from literature. 

The source code is in three parts. The file params . inc contains the grid 
dimensions while the file cavity . inp contains the input deck. The file cavity . f 
is the bulk of the code. To execute, compile cavity. f and run. This visual 
User’s Manual is derived from the source code via processing with the M£X2 r 
software. 

Following the Introduction are the two files which an' typically edited be- 
tween solutions, params. inc and cavity, inp. Then follows the driver routine 
and the subroutines, each as a separate section of the manual. 
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Params.inc 


The parameter file, params. inc, for cavity sets the number of points in i and 
j as imax and jmax. 

c parameter file for cavity params.inc 

parameter ( imax = 257, jmax = 257 ) 

parameter ( imxmed = (imax+l)/2, jmxmed = (jmax+l)/2 ) 
parameter ( imxcor = (imxmed+l)/2, jmxcor = ( jmxmed+1) /2 ) 

dimension dx( imax, jmax, 2), x(imax, jmax, 2), 

$ f(imax, jmax, 5), 

$ alam(imax, jmax, 2), sig(imax, jmax, 2) 

dimension fmed(imax, jmax, 5), fcor(imax, jmax, 5) 

dimension aij(imax, jmax), aipl(imax, jmax), aiml(imax, jmax), 

$ ajpl(imax, jmax), ajml(imax, jmax), b(imax, jmax), 

$ res (imax, jmax, 5, 3) 

common / vars / ntot , dx, x, f, alam, sig, fmed, fcor 
common / vecs / ai j , aipl, aiml , ajpl, ajml, b, res 


The dependent variables are stored in f(i,j,var) where var takes on values 
of 1-5 corresponding to (u, v , C, P, ij>). 
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Cavity, inp 

The input deck is cavity, inp. 

input deck to cavity 
$gr ids 
alength = 1 . 
height = 1. 
igrid = 0 
beta = 1.4 
$ 

$condit ions 
alammax = 1 . 
sigmax = 1. 
re = 2000. 

ubig = 1. 
igov = 2 
nitsgs = 1 
niterp = 1000 
convg = 5.e-5 
iskpout = 50 
img = 2 

$ 


c alength - cavity length 

c height - cavity height 

c igrid - 0 = uniform grid 

c 1 - beta clustering to edges 

c beta - clustering strength > 1 

c alammax - CFL-like number for inviscid terms 


c sigmax - CFL-like number for viscous terms 

c re - Reynolds number = ubig * alength / nu 

c ubig - upper boundary speed 
c igov - 1 = primitive variables formulation 
c 2 = stream f unction/vorticity formulation 

c nitsgs - number of sgs sub-sweeps 
c niterp - maximum number of iterations 
c convg - convergence tolerance 
c iskpout - convergence history output skip 
c img - 0 = fine mesh only 
c 1 = mesh sequencing 

c 2 = mesh sequencing + V-cycle multigrid 


cavity . inp 
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Main Driver 


The driver routine contains the logic for mesh sequencing and if multigrid 
is desired. Subroutines called: inputs, grid, dxdy, initial, timestep, 
solvr, restrict, be, prolong, multig, outsol. The subroutines follow 
in that order (the order in which they are called). 

This code was originally computer assignment 1 for Applied Science 778 at 
William & Mary, 1996. 

program cavity 

Load dimension, common, and parameter statements. 

include ’paran^in^ 
logical istop 

common / parms / alength, alammax, sigmax, re, height, ubig, 
$ igov, nitsgs, niterp, convg, iskpout, img, istop, 

$ igrid, beta 

common / parm2 / epustart, dxfin, dxmed, dxcor, iskpfin, 

$ iskpmed, iskpeor, residO, nout 

Default initial values. 

istop = .false, 
dxfin = 1. 
dxmed = . 5 

dxcor = .25 

iskpfin = 1 
iskpmed = 2 
iskpeor = 4 
ntot = 0 

Load input deck, cavity, inp. 
call inputs 

Set CPU timer at start of main routine, etime is timer for SUN and SGI. 
atimel, atime2 are dummy variables. 

epustart = etime ( atimel, atime2 ) 

Generate the fine grid. 

imxl = imax 

jmxl = jmax 

call grid ( igrid, imxl, jmxl, alength, height, beta ) 

Define dx and dy grid spacings on finest grid, 
call dxdy ( imxl, jmxl ) 


! Stopping flag 

! Fine grid 1:1 spacing, relative to fine 
! Medium grid 1:2 spacing, relative to fine 

! Coarse grid 1:4 spacing, relative to fine 

! Medium grid coarsed by 2 relative to fine 

! Coarse grid coarsed by 4 relative to fine 
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Initialize' fiowfield on finest grid. 


call initial ( imxl , jmxl, ubig, dxfin, iskpfin, f (1,1,1) ) 
Set dt time step size, and lambda, sigma coefficients. 

call timestep ( alammax, sigmax, re ) 

Run 1 timestep on fine mesh to get starting residual- res idO. 
residO = -1. 

call solvr ( 1, 0, imxl, jmxl, igov, dxfin, iskpfin, 

$ nitsgs, re, ubig, residO, istop, convg, nout , 

$ iskpout, epustart, f ( 1 , 1 , 1 ) ) 

Ski]) mesh sequencing if img=0. Perform mesh sequencing if img=l or 2. 

if ( img . eq. 0 ) go to 100 

Mesh Sequencing 

Restrict from fine to medium mesh. 

call restrict ( imxmed, jmxmed, f (1,1,1), fmed(l,l,l) ) 
Restrict from medium to coarse mesh. 

call restrict ( imxcor , jmxcor, fmed(l , 1 , 1) , fcor(l,l,l) ) 
Impose boundary conditions. 

call be ( igov, imxcor, jmxcor, ubig, re, iskpeor, dxcor, 

$ fcor(l,l,l) ) 

Iterate on coarse mesh for nitc iterations. Options for the number of iterations 
to perform ant 

ImaXc ■ Jmaj’r 
M>r = j 

N it c = ma x(Imax c . Jmax c ) 

Nit r = Irnax c + Jmax c 

nitc = imxcor * jmxcor / 2 
c nitc = max( imxcor, jmxcor ) 

c nitc = imxcor + jmxcor 

Limit the number of iterations on the coarse mesh to be no more than twice the 
maximum number of fine mesh iterations. 

nitc = min( nitc, niterp + niterp ) 

call solvr ( nitc, 0, imxcor, jmxcor, igov, dxcor, iskpeor, 

$ nitsgs, re, ubig, residO, istop, convg, nout, 

$ iskpout, epustart, f cor (1,1,1) ) 



Prolongate from coarse to medium mesh, impose boundary conditions, and relax 
on medium mesh. 

call prolong ( imxcor, jmxcor, fmed(l,l,l), fcor(l,l,l) ) 
call be ( igov, imxmed, jmxmed, ubig, re, iskpmed, dxmed, 

$ fmed(l ,1,1) ) 
nitc = imxmed * jmxmed / 2 
c nitc = max( imxmed, jmxmed ) 

c nitc = imxmed + jmxmed 

nitc = min( nitc, niterp + niterp ) 

call solvr ( nitc, 0, imxmed, jmxmed, igov, dxmed, iskpmed, 

$ nitsgs, re, ubig, residO, istop, convg, nout , 

$ iskpout, epustart, fmed( 1,1,1) ) 

Now prolongate from medium to fine mesh and impose boundary conditions. 
Ready to relax on fine mesh. 

call prolong ( imxmed, jmxmed, f (1,1,1), fmed(l,l,l) ) 
call be ( igov, imax, jmax, ubig, re, iskpfin, dxfin, 

$ f(l,l,l) ) 

Skip point for bypassing mesh sequencing. 

100 continue 

Relaxation on Fine Mesh 

Is multigrid desired? If so. img=2, and we go to subroutine multig to converge 
solution. Otherwise, we relax on the fine grid only. 

if ( img .eq. 2 ) then 
call multig 

else 

call solvr ( niterp, 0, imax, jmax, igov, dxfin, iskpfin, 

$ nitsgs, re, ubig, residO, istop, convg, nout, 

$ iskpout, epustart, f( 1,1,1) ) 

end if 

Output solution and end of program. 

call outsol ( imxl, jmxl, nout ) 

stop 

end 
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Inputs 

Inputs reads the input deck cavity, inp using the namelist construct. It is 
calk'd from cavity. 

subroutine inputs ! inputs 

common / parms / alength, alammax, sigmax, re, height, ubig, 
$ igov, nitsgs, niterp, convg, iskpout, img, istop, 

$ igrid, beta 

logical istop 

namelist / grids / alength, height, igrid, beta 
namelist / conditions / alammax, sigmax, re, ubig, igov, 

$ nitsgs, niterp, convg, iskpout, img 

open ( 9, file = ’ cavity . inp ’ , form = ’formatted’ ) 
read ( 9 , * ) 
read ( 9, grids ) 
read ( 9, conditions ) 

return 

end 
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Grid 


Grid is a 2-D grid generator for rectangular domains. The lower-left corner is 
set at the origin. Called by cavity. 

igrid = 0 — uniform spacing in x and y 

igrid = 1 Clustering to all walls, of strength 0 > 1. with stronger 
clustering as fi — > 1. 

subroutine grid ( igrid, imx, jrax, width, height, beta ) !..grid 

include * params . inc ’ 

Ax and Ay for uniform spacing. 

dxl =1. / float ( imx - 1 ) 
dyl =1. / float ( jmx - 1 ) 


0 + 1 and 0 — 1 . 


bpl = beta + 1 . 
bml = beta - 1 . 
d = .5 

Set x — 0 on the left-hand boundary and x — width on the right-hand boundary. 
x(i, j, 1) = ^-coordinate, x(i,j, 2) = //-coordinate. 

do j = 1, imx 
x(l, j,l) - 0. 
x (imx , j , 1) = width 
end do 

Set y = 0 on the lower boundary and y = height on the upper boundary. 

do i = 1 , imx 
x ( i , 1,2) = 0. 

x(i,jmx,2) = height 
end do 

Assign x values between left (/ = 1) and right (i — imx) sides, xl is the 
arclength fraction of the distance from the left to right for the current point. 

do i = 2, imx-1 

Either uniform spacing, 


Xi 


- Ax(i - 1) = 


i — 1 

Imau 1 


xl = dxl * float ( i - 1 ) 
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or beta clustering. 


T + R(Ii - 1) - (H + 1) 

Xl ~ T 

where. 

(ti+ i\ ( t=#) / — 1 

T = (2rf+l)(J? + l). R=(-—-) . A = - — 

if ( igrid . eq. 1 ) then 
an = xl 

rn = ( bpl / bml )**((an - d) / (1. - d) ) 
tn = ( 2 . * d + 1 . ) * ( rn + 1 . ) 
xl = ( tn + rn * bml - bpl ) / tn 
end if 

do j = 1, jmx 

x(i , j , 1) = x ( 1 , 1 , 1 ) + xl * width 
end do 
end do 

Fill in the y values between bottom (j = 1) and top ( j — jmx) using the same 
spacing variation as in the x direction. 

do j =2, jmx-1 

yl = dyl * float ( j - 1 ) 

if ( igrid .eq. 1 ) then 
an = yl 

rn = ( bpl / bml )**((an - d) / (1. - d)) 
tn = (2. * d + 1. ) * (rn+1. ) 
yl = ( tn + rn * bml - bpl ) / tn 
end if 

do i = 1 , imx 

x(i , j ,2) = x(l , 1 ,2) + yl * height 
end do 
end do 

return 

end 
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dxdy 

Computes Ax and Ay grid spacings on a cartesian-aligned rectangular grid. 
Second-order central differencing used for interior points. Second-order one- 
sided differencing used for boundary points. Called by cavity. 

subroutine dxdy ( imx, jmx ) ! dxdy 

include ’ params . inc 9 
At interior points Ax — - (;?•/+ 1 — j ) 
do i = 2 , imx - 1 

dxl = ( x(i+l,l,l) - x(i-l,l,l) ) * .5 
do j = 1, jmx 
dx (i , j , 1) = dxl 
end do 
end do 


At the left end Ax = ^(— 3.ri -I- 4x2 — # 3 ) 

dxl = ( -3. * x(l,l,l) + 4. * x(2,l,l) - x(3,l,l) ) * .5 
do j = 1 , jmx 
dx ( 1 , j , 1 ) = dxl 
end do 

At the tight end Ax — ^ (3x Ar — 1 “l - 3'inu— 2 ) 

dxl = ( 3 . * x(imx ,1,1) - 4.* x(imx-l,l,l) + x(imx-2,l,l) ) * .5 
do j = 1, jmx 

dx(imx , j ,1) = dxl 
end do 

Now do the same for Ay , interior, bottom, and top points, 
do j = 2, jmx - 1 

dyl = ( x(l , j+1 ,2) - x(l , j-1 ,2) ) * .5 
do i = 1 , imx 
dx (i , j ,2) = dyl 
end do 
end do 

Bottom. 

dyl = ( -3. * x (1 , 1 , 2) + 4. * x(l,2,2) - x(l,3,2) ) * ,5 
do i = 1 , imx 
dx(i,l,2) = dyl 
end do 


12 


Top. 


dyl = ( 3.* x(l,jmx,2) - 4.* x(l,jmx-l,2) + x(l, jmx-2,2) ) * .5 
do i = 1 , imx 

dx(i,jmx,2) = dyl 
end do 


return 

end 
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Initial 


Sets initial conditions for zero motion at t = 0 . At t = 0+ the upper surface 
impulsively started at speed U. Called by cavity. 

subroutine initial ( imx, jmx, u, dxfac, iskp, tf ) !... initial 

include } params . inc } 

tf is the array of solution variables to be initialized. 

dimension tf (imax , jmax , 5) 

First set everything to zero at t = 0 _ . 

do k = 1, 5 
do j = 1, jmx 
do i = 1 , imx 
tf (i, j ,k) = 0. 
end do 
end do 
end do 

Set the initial press un* to P rt j = 1 . 

do j = 1 , jmx 
do i = 1 , imx 
tf (i, j ,4) = 1. 
end do 
end do 

Now impulsively start the upper surface with speed U — ubig, set in cavity, inp. 
do i = 1 , imx 

tf(i,jmx,l) = u ! u-velocity 

tf(i,jmx,3) =2. * u * dxfac / dx(i,jmax,2) ! vorticity 

end do 

return 

end 
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Timestep 

timestep, railed from cavity, determines the At timestep size for either local 
or global timestepping. At is then incorporated into the parameters alam and 
sig. alamx and sigx are analogous to convective and diffusive CFL numbers, 
respectively, and are set in cavity, inp. 

subroutine timestep ( alamx, sigx, re ) !.... timestep 

include 5 params . inc ’ 

Initialize dtmin = min (A t) for global timestep option. 

dtmin = 2. * alamx * max( dx(2,2,l), dx(2,2,2) ) 

do j = 1, jmax 
do i = 1 , imax 

Compute the local convective At. 


At = 2A win . r min(A.r, Ay) 

dtl = 2. * alamx * min( dx(i,j,l), dx(i,j,2) ) 

Compute 1 the local diffusive At. 

At = R v a mtix min(A.r _) , Ay 2 ) 

dt2 = re * sigx * min( dx(i , j , 1) **2 , dx(i,j,2)**2 ) 

The local At is the minimum of these tiniest, eps. For global timestepping At is 
the minimum of all the local timesteps. 

dt = min( dtl, dt2 ) 
dtmin = min( dtmin, dt ) 

Now incorporate At into alam and sig coefficients. 

Ar 2 Ax ' u 2 Ay 

At _ At 

0i Re Ax' 2 ' ° XJ Re Ay 2 

do k = 1, 2 

alam(i,j,k) = dt * .5 / dx(i,j,k) 
sig(i , j ,k) = dt / re / dx (i , j ,k) **2 
end do 

end do 
end do 
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This section sets a global tirnestep for all points. If local timestepping is desired, 
comment out the following do loops. 

do j = 1, jmax 
do i = 1 , imax 
do k = 1, 2 

alam(i,j,k) = dtmin * .5 / dx(i,j,k) 
sig(i,j,k) = dtmin / re / dx(i,j,k)**2 
end do 
end do 
end do 

return 

end 
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Solvr 


Linear system solver for Af = b. Calk'd by cavity and multig. Matrix A is 
filled in subroutine loada and the right-hand side (RHS) vector b is assigned in 
subroutine loadb. Subroutines sgs. be. deep, and stopping are also called, and 
follow solvr in this printout. Af = b is solved iteratively in the time variable 
using symmetric Gauss-Seidel sweeps. 

mg = 0 working on top level of multigrid (fine grid) 
mg = 1 descending cycle of error smoothing 
mg = 2 ascending cycle of error smoothing 

subroutine solvr ( niterp, mg, imx, jmx, igov, dxfac, iskp, 

$ nitsgs, re, ubig, residO, istop, convg, nout , iskpout, 

$ epustart , tf ) ! . . 

include ? params.inc , 
logical istop 

dimension fold(imax, jmax, 5), tf (imax , jmax , 5) 

Evolve in time until convergence or niterp iterations (main loop). 

do n = 1, niterp 
ntot = ntot + 1 

Retain previous solution to check convergence rate. 

do k = 1, 5 
do j = 1, jmx 
do i = 1 , imx 

f old(i , j ,k) = tf (i , j ,k) 
end do 
end do 
end do 

Solve convection/diffusion equations. 

ft + Itfx + vfy - -i-V- / = -g 

First step is to load pent a- diagonal terms of A matrix in a, j . a j+ ] j . a i j , u, j+ ] , 
and (i{,j — i • 

call loada ( 1, imx, jmx, dxfac, iskp, mg, tf ( 1,1,1) ) 

If wo are solving the primitive* variables formulation. 

if ( igov .eq. 1 ) then 
Solve for u from. tt t + uu r 4- vu f/ = —P r + 
kk = 1 


solvr 
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Load RHS b vector if evolving on fine grid or the residual if this is a multi-grid 
error smoothing on a coarse grid. 

if (mg . eq. 0 ) then 

call loadb ( kk, imx, jmx, iskp, dxfac, “tf (1,1,1) ) 
else 

if ( mg . eq. 1 ) 11 = 1 

if (mg . eq. 2 ) 11 = 2 

do jj = 2, jmx-1 
do ii = 2, imx-1 

b(ii,jj) = res(ii , j j ,kk, 11) 
end do 
end do 
end if 

Now perform the symmetric Gauss-Seidel sweeps. 

call sgs ( kk, nitsgs, imx, jmx, tf (1,1,1) ) 

Save the residual. 

if (mg .eq. 0 ) 11 = 1 

if ( mg .eq. 1 ) 11 = 3 


do jj = 2, jmx- 1 
do ii = 2, imx-1 

res (ii , j j ,kk , 11) = b ( ii , j j ) 

$ - aij (ii , j j ) * tf (ii , j j ,kk) 

$ - aipl (ii , j j ) * tf (ii+1, j j ,kk) 

$ - aiml(ii,jj) * tf (ii-1 , j j ,kk) 

$ - ajpl (ii , j j ) * tf (ii, j j+l,kk) 

$ - ajml ( ii , j j ) * tf (ii , j j -1 ,kk) 


end do 
end do 

Now solve for v from, v t -f uv x + vv u = —P y + j^-S/ 2 v. 
kk=2 

if (mg .eq. 0 ) then 

call loadb ( kk, imx, jmx, iskp, dxfac, tf( 1,1,1) ) 
else 

if ( mg .eq. 1 ) 11 = 1 
if ( mg .eq. 2 ) 11 = 2 
do jj = 2, jmx-1 
do ii = 2, imx-1 

b(ii,jj) = res(ii , j j ,kk,ll) 
end do 
end do 
end if 
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call sgs ( kk, nitsgs, imx, jmx, tf(l,l,l) ) 


if (mg.eq. 0) 11= 1 
if ( mg .eq. 1 ) 11 = 3 
do jj = 2, jmx-1 
do ii = 2, imx-1 

res(ii, j j ,kk,ll) = b(ii,jj) 

$ - aij (ii , j j ) * tf (ii , j j ,kk) 

$ - aipl(ii,jj) * tf (ii+1 , j j ,kk) 

$ - aiml(ii,jj) * tf (ii-1 , j j ,kk) 

$ - ajpKii , j j) * tf (ii , j j+1 ,kk) 

$ - ajml(ii,jj) * tf (ii , j j -1 ,kk) 


end do 
end do 

If instead of primitive variables we are solving the stream-function and vortieity 
equations, we evolve ( from. 0 + + r Ci/ — 77- V 2 ( t . 

else if ( igov .eq. 2 ) then 
kk = 3 

if ( mg .eq. 0 ) then 

call loadb ( kk , imx, jmx, iskp, dxfac, tf( 1,1,1) ) 
else 

if ( mg .eq. 1 ) 11 = 1 
if ( mg .eq. 2 ) 11 = 2 
do jj = 2, jmx-1 
do ii = 2, imx-1 

b(ii , j j) = res (ii , j j , kk , 11) 
end do 
end do 
end if 

call sgs ( kk, nitsgs, imx, jmx, tf (1,1,1) ) 


if (mg .eq. 0 ) 11 = 1 
if (mg .eq. 1 ) 11 = 3 
do jj = 2, jmx-1 
do ii = 2, imx-1 

res (ii , j j ,kk , 11) = b ( ii , j j ) 

$ - aij (ii , j j ) * tf (ii , j j ,kk) 

$ - aipl(ii,jj) * tf (ii+1 , j j ,kk) 

$ - aiml(ii,jj) * tf (ii-1 , j j ,kk) 

$ - ajpl(ii,jj) * tf (ii , j j+1 ,kk) 

$ - ajml(ii,jj) * tf (ii , j j-1 ,kk) 


end do 
end do 

else 

write (6,*) * bad value for igov = igov 
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endif 


Having evolved the convective/diffusive equations, the Poisson equations remain 
to he solved. 

-V-7 = -y 

Load penta-diagonal terms for .4 matrix. 

call loada ( 2, imx, jmx, dxfac, iskp, mg, tf( 1,1,1) ) 

If using the primitive variables formulation we solve for P. 
if ( igov .eq. 1 ) then 

kk = 4 

if ( mg .eq. 0 ) then 

call loadb ( kk, imx, jmx, iskp, dxfac, tf ( 1,1,1) ) 
else 

if ( mg .eq. 1 ) 11 = 1 
if ( mg .eq. 2 ) 11 = 2 
do jj “ 2, jmx-1 
do ii = 2, imx-1 

b(ii,jj) = res(ii, j j ,kk, 11) 
end do 
end do 
end if 

call sgs ( kk, nitsgs, imx, jmx, tf (1,1,1) ) 


if (mg .eq. 0 ) 11 = 1 
if (mg .eq. 1 ) 11 = 3 
do jj = 2, jrax-1 
do ii = 2, imx-1 

res(ii, j j ,kk,ll) = b(ii.jj) 

$ - aij (ii , j j ) * tf (ii , j j ,kk) 

$ - aipl (ii, j j ) * tf (ii+1, jj ,kk) 

$ - aiml (ii , j j ) * tf (ii-1 , j j ,kk) 

$ - ajpl (ii , j j ) * tf (ii, j j+l,kk) 

$ - ajml (ii , j j ) * tf (ii , j j-1 ,kk) 


end do 
end do 

Otherwise solve for the stream-function 0 from. V* 2 0 = (. 

else 
kk = 5 

if (mg .eq. 0 ) then 

call loadb ( kk, imx, jmx, iskp, dxfac, tf ( 1,1,1) ) 
else 

if ( mg .eq. 1 ) 11 = 1 
if (mg .eq. 2 ) 11 = 2 
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do jj = 2, jmx-1 
do ii = 2, imx-1 

b(ii , j j) = res(ii , j j ,kk,ll) 
end do 
end do 
end if 

call sgs ( kk, nitsgs, imx, jmx, tf(l,l,l) ) 


if ( mg .eq. 0 ) 11 = 1 
if ( mg .eq. 1 ) 11 = 3 
do jj = 2, jmx-1 
do ii = 2, imx-1 

res(ii, j j ,kk,ll) = b(ii,jj) 

$ - aij (ii , j j ) * tf (ii , j j ,kk) 

$ - aipl (ii , j j ) * tf (ii+l,j j ,kk) 

$ - aiml(ii,jj) * tf (ii-1 , j j ,kk) 

$ - ajpl(ii, j j) * tf (ii, jj+l,kk) 

$ - ajml(ii,jj) * tf (ii, jj-l,kk) 


end do 
end do 
end if 

If we are performing an update on the fiat 1 mesh then update the boundary 
points. Also, for the stream-function/ vorticity formulation the velocities and 
pressure need to be decoded at the interior points. Otherwise, if this is an error 
smoothing iteration, skip this section. 

if (mg .eq. 0 ) then 

call be ( igov, imx, jmx, ubig, re, iskp, dxfac, 

$ tf (1,1,1) ) 

if ( igov .eq. 2 ) then 

Decode the velocities. 

do j = 2, jmx-1 

jl=(j-l)* iskp + 1 
do i =2, imx-1 

il=(i-l)* iskp + 1 

tf(i,j,l) = .5 * ( tf(i,j+l,5) - tf (i , j -1 , 5) ) 

$ / dx(il,jl,2) * dxfac 

tf (i , j ,2) = .5 * ( tf (i-1 , j ,5) - tf (i+1 , j ,5) ) 

$ / dx(il,jl,l) * dxfac 

end do 
end do 

Decode the pressure by iterating on the Poisson equation. 
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call deep ( imx, jmx, iskp, tf(l,l,l) ) 
end if 
end if 

Check convergence criteria to see if we stop. 

call stopping ( n, imx, jmx, igov, residO, istop, convg, 
$ resid, nout , iskpout, epustart, iskp, mg, 

$ fold(l,l,l), tf(l,l,l) ) 

Jump out of loop when converged. 

if ( istop .eq. .true. ) go to 100 

End of the main loop for time evolution. 

end do 

idebugwrite = 0 



if ( idebugwrite .eq. 

1 

) write 

(6,200) 

200 

format (/ ’reached maximum number of iterations without 1 , 


$ ’converging 5 /) 




100 

continue 





if ( idebugwrite .eq. 

1 

) write 

(6,201) n, iskp, resid 

201 

format ( 5 n = 5 , i5, 

> 

iskp = 5 

, i2, 5 resid = 5 , e!3 . 5 ) 


return 





end 
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Load A 


Form t ho peta-diagonal Left-Hand Side (LHS) implicit matrix. Stored in five 
diagonal vectors, Aij, A i±l j, and Ajj± j. loada is called from solvr. 

subroutine loada ( ieq, imx, jmx, dxfac, iskp, mg, tf ) 

include ’ params . inc ’ 
dimension tf (imax , jmax, 5) 


dl = dxfac 
d2 = dxfac**2 

For the convoctive/diffusive equations for a. i>, or 

Ajj = 1 H- 2(7 r 4- 2(7 y 

1 J — (7 x i \ r U . -4y , <7 y dl AjyC 

if ( ieq .eq. 1 ) then 
do j = 1, jmx 

jl=(j-l)* iskp + 1 
do i = 1 , imx 

il=(i-l)* iskp + 1 
if ( mg .eq. 0 ) then 
ul = tf (i, j ,1) 
vl = tf (i , j ,2) 
else 

ul = f (il, jl, 1) 
vl = f (il , j 1 ,2) 
end if 

aij (i , j)= 1. + 2 . *sig(il , j 1 , 1) *d2 + 2. *sig(il , j 1 ,2)*d2 
aipl(i,j) = -sig(il , j 1 , l)*d2 + alam(il , j 1 , 1) *dl * ul 
aiml (i , j ) = -sig(il , j 1 , 1) *d2 - alam(il , j 1 , 1) *dl * ul 
ajpl (i , j ) = -sig (il , j 1 , 2) *d2 + alam(il , j 1 , 2) *dl * vl 
ajml (i , j ) = -sig(il , j 1 ,2) *d2 - alam(il , j 1 , 2) *dl * vl 
end do 
end do 


For the Poisson equations for P or t/\ 


.4 


i'j 


2 2 
Ax 2 + A y- 


.4 


i±\J 


A./* 2 ' 


Aij±\ — - 


1 


else if ( ieq .eq. 2 ) then 
do j = 1, jmx 
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$ 


jl=(j-l)* iskp + 1 
do i = 1 , imx 

il = ( i - 1 ) * iskp + 1 
aij (i , j ) = 2 . * d2 / dx(il , j 1 , 1) **2 + 

2 . * d2 / dx(il, jl,2)**2 
aipl(i,j) = -d2 / dx(il , j 1 , 1) **2 
aiml (i , j ) = -d2 / dx(il , j 1 , 1) **2 
ajpl (i , j ) = -d2 / dx(il , j 1 ,2)**2 
ajml(i,j) = -d2 / dx(il , j 1 , 2) **2 
end do 
end do 

end if 

return 

end 
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Load b 


Assign the Right-Hand Side (RHS) vector b for advancing Af = b from time 
level u to level 7 / + 1. loadb is called by solvr. 

subroutine loadb ( mf , imx, jmx, iskp, dxfac, tf ) 
include ’ params.inc’ 
dimension tf (imax , jmax , 5) 


dl = dxfac 

d2 = dxfac * dxfac 


If solving for u 


if ( mf .eq. 1 ) then 
do j = 2, jmx - 1 

jl=(j-l)* iskp + 1 
do i = 2, imx - 1 

il = ( i - 1 ) * iskp + 1 

b(i,j) = tf (i , j ,mf ) - alam(il, jl,mf ) * dxfac * 
$ ( tf(i+l,j,4) - tf (i-1 , j ,4) ) 

end do 
end do 




If solving for i\ 


b = v - 




else if ( mf .eq. 2 ) then 
do j = 2, jmx - 1 

jl=(j-l)* iskp + 1 
do i = 2 , imx - 1 

il=(i-l)* iskp + 1 

b (i , j ) = tf (i , j ,mf ) - alam(il , j 1 ,mf ) * dxfac * 
$ ( tf (i , j+1 ,4) - tf (i , j-1 ,4) ) 

end do 
end do 


If solving for 

else if ( mf .eq. 3 ) then 
do j =2, jmx - 1 
do i = 2, imx - 1 
b(i, j) = tf (i,j ,mf ) 
end do 
end do 
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If solving for P, 


C (uj+i-Uj-iY 
''=~Ji + 5P — 


+ (Uj + 1 - «j-l )(«»+! -Vj-l) + (Vj + 1 - Ij-l) 


2A:r A?/ 


4Ay 2 


1 

7l 


g>+i -2C + C f --i Cj+i - 2C + C,_i 


Ax 2 


Ay 2 


where, 


c 

c 

c 

c 

c 

c 

c 

c 

c 


£1 iif— ) ^ Vj + 1 e j _ i 


2Ax 


2Ay 


else if ( mf . eq. 4 ) then 
do j = 2, jmx - 1 

jl=(j-l)* iskp + 1 
do i = 2, imx - 1 

il = ( i - 1 ) * iskp + 1 
b (i , j ) = .25 * d2*( (tf (i+1 , j , 1) - 
$ tf (i-1 , j , 1))**2 / dx(il , j 1 , 1)**2 

$ + 2 . * (tf (i , j+1 , 1) - tf(i, j-1,1)) * (tf (i+1 , j ,2) 

$ - tf (i-1 , j ,2)) / dx(il,jl,l) / dx(il , j 1 ,2) 

$ + (tf (i, j+1 ,2) - tf (i,j-l,2))**2 / 

$ dx(il, jl,2)**2 ) 

b(i , j ) = .5 * ( 

$ 

$ 

$ 

$ 

$ 

b(i,j) = .5 * 

$ 

$ 

end do 
end do 


( f (i, j + 1,1) - f(i, j-1,1) ) * 

( f ( i+1 , j ,2) - f (i-1 , j ,2) ) - 
( f (i+1, j ,1) - f (i-1, j ,1) ) * 

( f(i,j+l,2) - f (i , j-1 ,2) ) ) 

/ dx ( i , j ,1) / dx(i , j ,2) 

( f (i+1 , j , 1) - f (i-1 , j , 1) ) **2 / 
dx(i , j , 1) **2 + .5 * ( f(i, j+1,1) - f(i, j-1,1) ) * 

( f(i+l,j,2) - f(i-l,j,2) ) / dx(i,j,l) / dx(i,j,2) 


If solving for psL b = — C- 

else if ( mf .eq. 5 ) then 
do j = 2, jmx - 1 
do i = 2 , imx - 1 

b(i , j ) = - tf (i, j ,3) 
end do 
end do 


else 

write (6,*) ’ mf not . le. 5 in loadb, mf = ’ , mf 
end if 
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c * diagnostic * 

c do i = 2, imx-1 

c do j = 2, jmx-1 

c write (6,*) 7 b^jbCijj), 7 for i 

c end do 

c end do 

return 

end 
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Sgs 

Symmetric GauB-Seidel iteration module. Iterates on Af = fc, sweeping first 
on increasing i nested inside increasing j , then sweeps again on decreasing j 
nested inside decreasing i. niter symmetric sweeps are performed. The vari- 
able/equation to be updated is designated by mf . sgs is called by solvr and 
deep. 


subroutine sgs ( mf , niter, imx, jmx, tf ) 

include J params . inc } 
dimension tf (imax , jmax , 5) 

imxl = imx - 1 
jmxl = jmx - 1 

How many symmetric sweeps to perform — niter . 
do k = 1, niter 

Forward sweep / and j increasing. Since .4 is penta-diagonal, the update 

formula is, 

fi'j ~ — (^ij — — (li-ijfj-ij — fi j-j-i — 

do j = 2, jmxl 
JPl = J + 1 
jml = j - 1 
do i =2, imxl 
ipl = i + 1 
iml = i - 1 

tf (i , j ,mf ) = ( b(i, j) - aipl(i, j) * tf (ipl , j ,mf ) - 
$ aiml(i, j) * tf(iml,j,mf) - ajpl(i,j) * tf(i,jpl,mf) 

$ - ajml(i, j) * tf (i , jml ,mf ) ) / aij(i,j) 

end do 
end do 

Backward sweep — i and j decreasing. 

do i = imxl, 2, -1 
ipl = i + 1 
iml = i - 1 
do j = jmxl, 2, -1 

jpl = j + 1 
jml = j - 1 

tf (i , j ,mf ) = ( b(i, j) - aipl(i, j) * tf(ipl,j,mf) - 
$ aiml(i, j) * tf(iml,j,mf) - ajpl(i,j) * tf(i,jpl,mf) 

$ - ajml (i , j ) * tf (i , jml ,mf ) ) / aij(i,j) 
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end do 
end do 

end do 

return 

end 
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Boundary Conditions 

Explicitly update the boundary points. All boundaries are solid, no-slip walls. 
Upper boundary translates at speed ubig. be is called by solvr. 

subroutine be ( ieq, imx, jmx, ubig, re, iskp, dxfac, tf ) 

include 5 params . inc , 
dimension tf (imax , jmax , 5) 

d2 = dxfac * dxfac 
i = 1, left wall 
i = 1 

il=(i-l)* iskp + 1 
do j =2, jmx-1 

jl = ( j - 1 ) * iskp + 1 

c if ( ieq . eq. 1 ) then ! primitive variables formulation 

a = v = 0 


tf (i,j ,1) = 0. 
tf (i, j ,2) = 0. 

Several options for pressure* are: 


P wall — Pwail+\ 


C 

c 

c 

c 

c 

c 

c 

c 


4 1 

Pwntl = ’qPwuU+I ~ trait +2 


f*wait — ^ ( i H*tvnll -\- 1 P-wull+2 4 “ 


4-Pu?rt//+ 1 “ 2P, 


R r Ax 


“ ,fW/ +2 9 p 2 \ 

+ *^waU+ 1 J 


ry ry . 1 ^U’a//+2 . ry'2 

^ wall — *wall + 1 + RAx ^ *watl+l 


tf (i, j ,4) = tf (i+1, j ,4) 

f (i , j ,4) * 4. / 3. * f (i+1 , j ,4) - f (i+2 , j ,4) / 3. 
f (i , j ,4) = ( 4. * f (i+1 , j ,4) - f (i+2 , j ,4) + 

$ ( 4. * f (i+1 , j , 1) - 2. * f (i+2, j , 1) ) 

$ / re / dx ( i , j , 1 ) + 2. * f (i+1 , j , 1 ) **2 ) / 3. 

f(i,j,4) = f (i+1 , j ,4) + ( 2. * f (i+1 , j , 1) - f (i+2, j , 1) ) 
$ / re / dx(i,j,l) + f (i+1 , j , 1) **2 


else 


vorticity/stream function formulation 


Options for i/t and ( are: 


V’wall — 0 
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c 

c 


wnl I — b'wall - fl 

_ 4 1 

C'i rail — T^y’wall+i ^ 4 <<j //-)-2 


(.wall 


2li'waU -\- 1 

A:r 2 


Clt'ttll = "7 — ~ (V’luall ~ 2fc’, , ,,//+ | + t’„, „»+•_>) 
A./— 


tf (i,j ,5) = 0. 

f (i , j ,5) = 4. / 3. * f (i+1, j ,5) - f (i+2, j ,5) / 3. 
f (i, j ,5) = f (i+l,j ,5) 


tf (i, j ,3) = 2. * d2 * tf (i+1 , j ,5) / dx(il , j 1 , 1) **2 
c f (i , j ,3) = ( f (i , j ,5) - 2. * f (i+1 , j ,5) + f(i+2,j,5) ) / 

c $ dx(i , j , 1)**2 

c end if 

end do 

i = imx, right wall 


i = imx 

il=(i-l)* iskp + 1 
do j = 2, jmx-1 

jl=(j-l)* iskp + 1 
c if ( ieq . eq. 1 ) then 

tf (i , j , 1) = 0. 
tf (i, j ,2) = 0. 

c f (i , j ,4) = 4. / 3. * f(i-l,j,4) - f (i-2 , j ,4) / 3. 

c f (i , j ,4) = f (i-1 , j ,4) - ( 2. * f (i-1 , j , 1) - f(i-2,j,l) ) 

c $ / re / dx ( i , j , 1 ) + f (i-1 , j , 1)**2 

tf(i,j,4) = tf(i-l,j,4) 
c else 

c f (i, j ,5) = 4. / 3. * f (i-1 , j ,5) - f(i-2,j,5) / 3. 

c f (i, j ,5) = f (i-1, j ,5) 

tf (i , j ,5) = 0 . 

c f (i, j ,3) = ( f (i , j ,5) - 2. * f (i-1 , j ,5) + f(i-2,j,5) ) / 

c $ dx(i , j , 1) **2 

tf (i, j ,3) =2. * d2 * tf (i-1 , j ,5) / dx(il , j 1 , 1)**2 
c end if 

end do 

j = 1, bottom 

j = i 

jl=(j-l)* iskp + 1 
do i = 2 , imx-1 

il = ( i - 1 ) * iskp + 1 
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c 


if ( ieq .eq. 1 ) then 
tf (i,j ,1) = 0. 
tf (i, j ,2) = 0. 

c f(i,j,4) = 4. / 3. * f (i,j+l,4) - f (i , j+2,4) / 3. 

c f (i , j ,4) = f (i , j+1 ,4) + ( 2. * f (i , j+1 ,2) - f(i,j+2,2) ) 

c $ / re / dx(i,j,2) + f (i , j+1 ,2) **2 

tf (i,j ,4) = tf (i, j+1, 4) 

c else 

c f (i » j , 5) = 4. / 3. * f (i, j+1 ,5) - f(i,j+2,5) / 3. 

c f (i.j ,5) = f (i,j+l,5) 

tf (i, j ,5) = 0. 

c f (i , j , 3) = ( f (i , j ,5) - 2. * f (i, j + 1,5) + f(i,j+2,5) ) / 

c $ dx(i,j,2)**2 

tf(i,j,3) =2. * d2 * tf(i,j+l,5) / dx(il, jl,2)**2 
c end if 

end do 

j = jinx, top 


J = jmx 

jl = (j-l)* iskp + 1 
do i =2, imx-1 

il=(i-l)* iskp + 1 
c if ( ieq .eq. 1 ) then 

Same options for v and P. u = U (ubig). 


V*watt — ^ Hawaii — {^V'wall — 1 


- Vhvait- > + 2 AyU) 


C wall — a •} (V J wall— 1 4" C wall 

Ay- 



W 


wail 


1 4" Wri// — 2 


tf (i, j ,1) = ubig 
tf (i, j ,2) = 0. 

c f(i,j,4) = 4. / 3. * f (i , j-1 ,4) - f(i,j-2,4) / 3. 

c f(i,j,4) = f (i , j-1 ,4) - ( 2. * f (i , j-1 ,2) - f(i,j-2,2) ) 

c $ / re / dx(i , j ,2) + f (i , j-1 ,2) **2 

tf (i, j ,4) = tf (i, j-1, 4) 
c else 

c f (i.j .5) = ( 4. * f (i , j-1 ,5) - f (i , j-2,5) + 

c $ 2. * ubig * dx(i,j,2) ) / 3. 

tf (i , j ,5) = 0. 

c f (i , j .3) = ( f (i.j >5) - 2. * f (i , j-1 ,5) + f(i, j-2,5) ) / 

c $ dx(i,j,2)**2 

tf (i, j ,3) = 2. * d2 * ( tf (i , j-1 ,5) + 

$ ubig * dx(il,jl,2) / dxfac ) / dx(il , j 1 ,2) **2 

c end if 

end do 
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just set the corner points to a neighboring point 

do k = 1 , 5 

tf (1 , 1 ,k) = tf (2, 1 ,k) 
tf(imx,l,k) = tf (imx-1 , 1 ,k) 
tf(l,jmx,k) = tf(2,jmx,k) 
tf (imx, jmx ,k) = tf (imx- 1 , jmx , k) 
end do 

return 

end 
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Deep 

Decode the pressure from the stream-function and vortidty. Relies on the fact 
that the velocities have already been decoded in solvr. The Poisson equation 
for the steady-state pressure, 

^ P — 2 Vy HyV,v ) 

is relaxed using symmetric Gaufi-Seidel sweeps. 

Called from cavity. Calls sgs. 

subroutine deep ( imx, jmx, iskp, tf ) 
include ’ params.inc* 
dimension tf( imax, jmax, 5 ) 

Build pentadiagonal LHS and forcing function RHS at interior points. Multiply 
through by grid spacing for convenience, giving. 

Ax 2 Ay 2 S7P = 2 Ax 2 Ay 2 ( u x v y — u y v x ) 

do i = 2 , imx - 1 

il = ( i - 1 ) * iskp + 1 
do j = 2, jmx - 1 

jl=(j-l)* iskp + 1 

Construct LHS with two-dimensional, second-order central differences. 

Ax 2 Ay 2 V 2 Pj j -> Ar(P,+ i + Pi- 1 ) + Ax 2 (Pj+i + Pj^ ) - 2(Ax 2 + Ay 2 )P UJ 

aij (i , j) = -2. * ( dx (il , j 1 , 1) **2 + dx(il , j 1 ,2) **2 ) 

aipl (i , j ) = dx(il,jl,2)**2 

aiml(i,j) = dx (il , j 1 ,2) **2 

ajpl(i,j) = dx(il, jl,l)**2 

ajml (i , j ) = dx(il, jl,l)**2 

Construct the RHS vector using second-order central differences. 


2 Ax' 2 Ay 2 ( u x v y - u y v r ) 


“ Ui-\)(v J+l - vj _ ] ) - («j+i - - c/_i)] 


ux = f (i+1, j , 1) - f (i-1 , j , 1) 
uy = f (i , j+1 , 1) - f(i,j-l,l) 
vx = f (i+1 , j ,2) - f (i-1 , j ,2) 
vy = f (i, j+1,2) - f (i , j-1 ,2) 
b(i , j) = .5 * dx(il,jl,l) * dx(il , j 1 ,2) * 
$ ( ux * vy - uy * vx ) 

end do 
end do 


34 



Finally iteratively solve for the pressure using Gaufi-Seulel sweeps. 

call sgs ( 4, 1, imx, jmx, tf(l,l,l) ) 

return 

end 
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Stopping 

Check to see if convergence criteria lias been met. Also outputs convergence 
and timing history. Called by solvr. 

subroutine stopping ( n, imx, jmx, igov, residO, istop, 

$ convg, resid, nout , iskpout, cpustart, iskp, mg, 

$ fold, tf ) 


include ’params.inc* 

dimension f old (imax , jmax , 5) , tf (imax , jmax , 5) 
logical istop 

Compute L‘2 residual of all variables at all points. 


resid = 0. 
do j = 1, jmx 
do i = 1 , imx 

c primitive variables 

if ( igov .eq. 1 ) then 

resid = resid + ( tf(i,j,l) 
resid = resid + ( tf(i,j,2) 
resid = resid + ( tf(i,j,4) 
else 

c stream-f unct ion/vorticity 

resid = resid + ( tf(i,j,3) 
resid = resid + ( tf(i,j,5) 
end if 
end do 
end do 

resid = sqrt( resid ) 

Initialize starting residual if first, iteration. 


fold(i,j,l) ) **2 
f old(i , j ,2) ) **2 
fold(i,j,4) )**2 


f old(i , j ,3) ) **2 
f old(i , j ,5) ) **2 


if ( residO .It. 0. ) then 
residO = resid 

write (45) 1, etime (al , a2) -cpustart , 1. 
write (6,*) J residO = ' , residO 
nout = 1 
end if 


Normalize the residual by the starting residual. 

resid = resid / residO 

Check for convergence if on fine mesh. 

if ( iskp .eq. 1 .and. abs (resid) .It. convg ) then 
istop = .true. 

write (6,*) 3 *** Solution Converged ***> 
end if 
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Save convergence history if on fine mesh. 

if ( ( mod( ntot, iskpout ) .eq. 0 
$ .or. istop .eq. .true. ) .and. mg .eq. 0 )then 

cpul = etime( atimel, atime2 ) - cpustart 
write(45) ntot, cpul, resid 
write (6,110) ntot, iskp, cpul, resid 
110 format ( 'iter = i5, ’ , iskp = ’ , i2,’, cput = ’ , 

$ f7.1, 's,', ' resid = ' , el2.6 ) 

nout = nout + 1 
end if 

return 

end 
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Restrict 


Restrict a fine-mesh vector to a coarse mesh by simple injection of the overlap- 
ping points, vbig is the vector to be restricted, defined on the fine grid, vsmal 
is the resultant vector, defined on the coarse grid, ix and jx are dimensions of 
vsmal. Called from cavity and multig. 

subroutine restrict ( ix, jx, vbig, vsmal ) 
include 5 params . inc y 

dimension vbig(imax , jmax , 5) , vsmal(imax, jmax ,5) 

do k = 1, 5 
dO j « 1, jX 

j2 = 2 * j - 1 
do i = 1 , ix 
i2 = 2 * i - 1 

vsmal (i , j ,k) = vbig(i2, j2,k) 
end do 
end do 
end do 

return 

end 
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Prolong 

A linear-averaging prolongation of a vector from a coarse grid to a fine grid, 
vsmal is the vector to prolongate, defined on the coarse grid, vbig is the result 
of the prolongation, defined on the fine grid, ix and jx are dimensions of vsmal. 
Called by cavity and multig. 

subroutine prolong ( ix, jx, vbig, vsmal ) 
include > params.inc’ 

dimension vbig(imax , jmax , 5) , vsmal(imax, jmax , 5) 
do k = 1, 5 

Simple injection for the overlapping points. 

do i = 1, ix 
j2 = 2 * j - 1 
do i = 1 , ix 
i2 = 2 * i - 1 

vbig(i2 , j2 ,k) = vsmal(i,j,k) 
end do 
end do 

Two-point averaging in i for the new fine-grid points along a j = const line of 
tin 1 roarse grid. 

do j = 1, jx 
j2 = 2 * j - 1 
do i = 1 , ix-1 
i2 = 2 * i 

vbig(i2, j2,k) = . 5 * ( vsmal (i , j ,k) + vsmal(i+l , j ,k) ) 
end do 
end do 

Two-point averaging in j for the new fine-grid points along a i = const line of 
the coarse grid. 


do i = 1, ix 

i2 = 2 * i - 1 
do j = 1, jx-1 
j2 = 2 * j 

vbig(i2, j2,k) = .5 * ( vsmal (i , j ,k) + vsmal (i , j+1 ,k) ) 
end do 
end do 
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Four-point averaging for remaining fine-grid points, corresponding to i ± j ± ^ 
on the coarse grid. 

do j = 1, jx - 1 
j2 = 2 * j 
do i = 1, ix - 1 
i2 = 2 * i 

vbig(i2, j2,k) = . 25 * ( vsmal(i,j,k) + vsmal(i+l »j,k) + 

$ vsmal(i, j+l,k) + vsmal (i+1 , j+1 ,k) ) 

end do 
end do 

end do 

return 

end 
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Multig 

V-cyde multigrid driver. Called from cavity. Calls solvr. restrict, and 
prolong. 

subroutine multig 

include ’params . inc ’ 
logical istop 

common / parms / alength, alammax, sigmax, re, height, ubig, 
$ igov, nitsgs, niterp, convg, iskpout, img, istop, 

$ igrid, beta 

common / parm2 / cpustart , dxfin, dxmed, dxcor, iskpfin, 

$ iskpmed, iskpcor, residO, nout 

Begin main loop of Y-cvcle. P(*rform niterp cycles. 

do nn = 1 , niterp 

Relax solution on fine-mesh. First, argument, is the number of smoothing itera- 
tions to perform. 

call solvr ( 3, 0, imax, jmax, igov, dxfin, iskpfin, 

$ nitsgs, re, ubig, residO, istop, convg, nout, 

$ iskpout, cpustart, f( 1,1,1) ) 

Are we converged on the fine mesh? 

if ( istop .eq. .true. ) go to 200 

Restrict residual from fine mesh to medium mesh. 

call restrict (imxmed,jmxmed, res (1,1,1, 1) ,res(l , 1 , 1 ,2) ) 
do k = 1, 5 

do j = 1, jmxmed 
do i = 1 , imxmed 

res(i,j,k,l) = res(i,j,k,2) 

fmed(i,j,k) =0. ! zero out error on medium mesh 

end do 
end do 
end do 

Smooth error on medium mesh. Again, the first argument is the number of 
smoothing iterations to perform. 

call solvr (2, 1, imxmed, jmxmed, igov, dxmed, iskpmed, 

$ nitsgs, re, ubig, residO, istop, convg, nout, 

$ iskpout, cpustart, fmed( 1,1,1) ) 

Restrict residual from medium to coarse mesh. 
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call restrict ( imxcor , j mxcor , res ( 1 , 1 , 1 , 3) , res ( 1 , 1 , 1 , 1 ) ) 
do k = 1 , 5 

do j = 1, jmxcor 
do i = 1 , imxcor 
f cor (i , j ,k) = 0 . 
end do 
end do 
end do 

Smooth error on coarse mesh. 

call solvr (1, 1, imxcor, jmxcor, igov, dxcor, iskpcor, 
$ nitsgs, re, ubig, residO, istop, convg, nout, 

$ iskpout, cpustart, fcor(l,l,l) ) 

Prolong coarse-mesh error to medium mesh. 

call prolong (imxcor, jmxcor, res (1 , 1 , 1 ,3) , fcor(l,l,l)) 

Correct medium-mesh error with contribution from coarse mesh. 

do k = 1, 5 

do j =2, jmxmed-1 
do i =2, imxmed- 1 

fmed(i,j,k) = fmed(i,j,k) + res(i,j,k,3) 
end do 
end do 
end do 

Smooth medium-mesh error again. 

call solvr (2, 2, imxmed, jmxmed, igov, dxmed, iskpmed, 
$ nitsgs, re, ubig, residO, istop, convg, nout, 

$ iskpout, cpustart, fmed( 1,1,1) ) 

Prolong medium-mesh error to fine mesh. 

call prolong (imxmed, jmxmed, res (1 , 1 , 1 , 1) , fmed(l , 1 , 1)) 

Update 1 fine- mesh solution with correction to residual from medium mesh. 

do k = 1, 5 

do j =2, jmax-1 
do i =2, imax-1 

f ( ± , j , k) = f(i,j,k) + res (i , j ,k , 1) 
end do 
end do 
end do 

Smooth fine- mesh solution again. 
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call solvr ( 3, 0, imax, jmax, igov, dxfin, iskpfin, 
$ nitsgs, re, ubig, residO, istop, convg, nout , 

$ iskpout, cpustart, f( 1,1,1) ) 

Now loop back to repeat V-cycie. 

end do 

200 continue 

return 

end 
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Outsol 


Solution-output routine* called from cavity. Solution is written to cavity_tec . dat 
in Tecplot[2] format for variables: x, y, Ax. Ay * ?r, v, P, u x + v u . 

Also outputs plot file of convergence history to cavity_convg.dat. 

subroutine outsol ( imx, jmx, nout ) 

include ’ params * inc J 
dimension deriv (imax , jmax ,2) 

Compute derivatives of velocities to verify continuity equation. Central differ- 
ence interior points, one-sided differences on boundaries. Should get u x +v y = 0. 

Ux 


do j = 1, jmx 

deriv(l , j , 1) = ( f(2,j,l) - f(l,j,l) ) / dx(l,j,l) 
do i = 2, imx-1 

deriv(i , j , 1) = .5 * ( f(i+l,j,l) - f (i-1 , j , 1) ) 

$ / dx(i,j ,1) 

end do 

deriv (imx, j , 1) = ( f(imx,j,l) - f (imx-1 , j , 1) ) 

$ / dx(imx, j , 1) 

end do 


do i = 1 , imx 

deriv(i ,1,2) = ( f(i,2,2) - f(i,l,2) ) / dx(i,l,2) 
do j = 2, jmx-1 

deriv(i, j ,2) = . 5 * ( f(i,j+l,2) - f(i,j-l,2) ) 

$ / dx(i, j ,2) 

end do 

deriv(i , jmx ,2) = ( f(i,jmx,2) - f (i, jmx-1, 2) ) 

$ / dx(i,jmx,2) 

end do 

Tecplot-format grid and solution file. 

open ( 15, file = ’cavity_tec.dat’, form = ’formatted’ ) 
write ( 15, 100 ) imx, jmx 
100 format ( ’TITLE = "flowsol"’ / 

$ ’VARIABLES = "X", "Y" , "dx", "dy"’/ 

$ ’"u", "v", "zeta",’ / 

$ ”’P M , "psi" , "cont"’ / 

$ ’ZONE 1=’, i3 , ’, J=’, i3, ’, F=BL0CK’ ) 

write ( 15, * ) ((( x(i,j,k), i=l,imx), j=l,jmx), k=l,2) 
write ( 15, * ) ((( dx(i,j,k), i=l,imx), j=l,jmx), k=l,2) 
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write ( 15, * ) ((( f(i,j,k), i=l,imx), j = l,jmx), k=l,5) 
write ( 15, * ) ((( deriv(i,j,l) + deriv(i,j,2) ), 

$ i=l,imx), j=l,jmx) 

Tecplot-format. convergence history. 

open ( 16, file = ’ cavity_convg . dat ’ , form = ’formatted' ) 
write ( 16, 101 ) nout 
101 format ( ’TITLE = "Convergence"’ / 

$ ’VARIABLES = "iter", "CPU", "R-L_2"’ / 

$ ’ZONE 1=’, i5 , J=l, F=point ’ ) 

rewind (45) 
do n = 1 , nout 

read (45) i, cpu, resid 
write (16,*) i, cpu, resid 
end do 

return 

end 
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Results 

Incompressible Navier-Stokes solutions were computed inside a unit cavity, with 
the upper surface translating to the right. Three Reynolds numbers ( R e ) were 
considered, 100, 400, and 2000. In all cases uniform meshes were employed, 
with the same spacing in both x and y directions. 

Convergence times for six cases, corresponding to a coarse and fine mesh 
convergence study at each of the three Reynolds numbers, are tabulated in 
Table 1. Solution times are seen to increase with increasing R e , and to increase 
by an order of magnitude with a quadrupling of the grid density. Solution times 
reported are for an SGI 195 MHz R10000 CPU. Each case was started from 
a stationary initial condition, using three-level mesh sequencing and multigrid 
smoothing to accelerate convergence*. 


R t 

,J 

CPU sec. 

L-2 

Resolution 

100 

51 

13 

HP 5 

±0.01 

100 

101 

121 

10 -6 

±0.005 

400 

101 

168 

10~ 6 

±0.005 

400 

126 

224 

1()- 6 

±0.004 

2000 

125 

707 

2 x 10 _r ’ 

±0.004 

2000 

249 

4039 

5 x 1(T 5 

±0.002 




Ghia 

Schreiber 

Present (coarse) 

Present (fine) 

J’o 

0.617 

0.G17 

0.62 

0.615 

yo 

0.734 

0.742 

0.74 

0.740 

Vo 

-0.103 

-0.103 

-0.102 

-0.103 

Co 

3.17 

3.18 

3.18 

3.20 

.ri 

0.031 

0.033 

0.03 

0.035 

Vi 

0.039 

0.025 

0.03 

0.035 

V’i (xl() e ) 

1.74 

2.05 

3.22 

2.09 

Cl 

-0.0255 

-0.0080 

-0.0231 

-0.0153 

Lx 

0.078 

— 

0.09 

0.085 

H x 

0.078 


0.09 

0.085 

x> 

0.945 

0.942 

0.94 

0.940 

v 2 

0.063 

0.050 

0.06 

0.060 

■t/b (xlO 5 ) 

1.25 

1.32 

1.59 

1.33 

C‘_> 

-0.0331 

-0.0255 

-0.0356 

-0.0359 

L-2 

0.133 


0.13 

0.135 

H> 

0.148 


0.15 

0.155 


Table 2: Vortex locations, stream-function, and vorticity for B f = 100. 



Figure 2: Vortex center arid separation point naming conventions. 
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At R t = 2000, six vorticies are predicted, seen in Figure 1. In addition to 
the primary vortex, there are three secondary vortices in the lower corners and 
upper-left side. The two remaining tertiary vortices are located in the extreme 
lower corners. Figure 2 is a cartoon sketch identifying the nomenclature for 
describing the vortices and their strengths. 

Numerical data from the six cases in Table 1 in the form of vortex locations, 
stream-function and vorticity strengths, and separation points are compared 
with the results of Ghia et a/[3] and Schreiber and Keller [4] in Tables 2 and 3. 
Excellent agreement is seen for R f = 100, Table 2, between the present method 
and the other sources. At this Reynolds number only three vortices are present. 

Table 3 contains the data for R e = 400. Schreiber reports data for only three 
vortices at this Reynolds number, though the present method resolved four and 
Ghia reports five vortices. Again, excellent agreement is seen between the data 
sets. Table 4 lists the results of the grid-convergence study at R e = 2000, where 
all six vortices appear. 

Summary 

An annotated source code for the solution of incompressible viscous flow in 
a two-dimensional cavity is presented. The Navier-Stokes equations are cast 
in a stream-function/ vorticity formulation with second-order spatial accuracy. 
Convergence to a steady state is accelerated via multigrid relaxation. 

Solutions are obtained for Reynolds numbers of 100, 400. and 2000. Solution 
accuracy is checked versus previously published numerical results. Excellent 
agreement is seen in the location of vortex centers, separation points, vorticity 
strengths, and stream-function values amongst the data sets. 

The present code is seen to be an accurate and efficient means for solving 
the incompressible viscous driven cavity problem. 
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Ghia 

Schreiber 

Present (coarse) 

Present (fine) 

r 0 

0.555 

0.557 

0.555 

0.556 

.Vo 

0.606 

0.607 

0.610 

0.608 

V'o 

-0.114 

-0.113 

- 0.112 

-0.113 

Co 

2.29 

2.28 

2.27 

2.29 


0.051 

0.050 

0.050 

0.052 

Vi 

0.047 

0.043 

0.045 

0.048 

fi (xlO 5 ) 

1.42 

1.45 

1.47 

1.49 

Cl 

-0.0570 

-0.0471 

-0.0600 

-0.0534 

L , 

0.127 


0.125 

0.128 

H i 

0.108 

- 

0.100 

0.108 

J'-J 

0.891 

0.886 

0.885 

0.884 

Vt? 

0.125 

0.114 

0.125 

0.124 

(.•■_> { X 10 4 ) 

6.62 

6.44 

6.43 

6.45 

Cj 

-0.434 

-0.394 

-0.449 

-0.455 

L> 

0.262 


0.265 

0.264 

H-2 

0.320 


0.320 

0.324 

J'4 

0.004 


... 


V-l 

0.004 

- 



i/\i (xlO 10 ) 

-7.67 




Ci 

0.0009 



-■ 


0.004 





0.004 


-- 



0.992 


0.990 

0.992 

Hr, 

0.008 


0.010 

0.008 

V’r, (xlO K ) 

- 1.86 


- 10.0 

-7.76 

Cs 

0.0044 


0.0092 

0.0052 

Lr, 

0.016 

— 

0.015 

0.016 

H 5 

0.016 

- 

0.015 

0.016 


Table 3: Vortex locations, stream-function, and vorticity for R, = 400. 
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Present (coarse) 

Present (fine) 

Xo 

0.524 

0.522 

Vo 

0.548 

0.556 

0o 

-0.120 

-0.134 

Co 

2.17 

2.57 

^1 

0.089 

0.087 

m 

0.097 

0.096 

V»i (xlO 5 ) 

66.6 

51.5 

Cl 

-0.679 

-0.615 

L x 

0.274 

0.270 

Hi 

0.214 

0.207 

X 2 

0.843 

0.833 

y 2 

0.099 

0.102 

V>-> (xl()') 

24.3 

24.7 

<2 

-1.67 

-1.65 

L, 

0.338 

0.360 

H-2 

0.395 

0.403 


0.033 

0.032 

V :i 

0.884 

0.884 

(XlO 1 ) 

1.47 

1.50 

Cl 

-0.806 

-0.943 

m 

0.073 

0.073 

Hi 

0.226 

0.224 

x A 

0.008 

0.006 

!M 

0.008 

0.006 

V>4 (XlO 8 ) 

-10.8 

-2.83 

c< 

0.0153 

0.0056 

L.\ 

0.012 

0.012 

Hi 

0.012 

0.012 

Xfi 

0.992 

0.991 

y 5 

0.008 

0.009 

V 5 (XlO 8 ) 

-25.6 

-11.5 

Cr> (xlO 3 ) 

9.89 

7.37 

L 5 

0.024 

0.022 

H 5 

0.024 

0.022 


Table 4: Vortex locations, stream-function, and vorticity for Ti, - 2000. 
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